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ABSTRACT 

The effects of a non-gradient flux term originating from the motion of convective elements with entropy 
perturbations of either sign are investigated and incorporated into a modified version of stellar mixing length 
theory (MLT). Such a term, first studied by Deardorff in the meteorological context, might represent the effects 
of cold intense downdrafts caused by the rapid cooling in the granulation layer at the top of the convection 
zone of late-type stars. These intense downdrafts were first seen in the strongly stratified simulations of Stein 
& Nordlund in the late 1980s. These downdrafts transport heat nonlocally, a phenomenon referred to as entropy 
rain. Moreover, the Deardorff term can cause upward enthalpy transport even in a weakly Schwarzschild-stably 
stratified layer. In that case, no giant cell convection would be excited. This is interesting in view of recent 
observations, which could be explained if the dominant flow structures were of small scale even at larger depths. 
To study this possibility, three distinct flow structures are examined: one in which convective structures have 
similar size and mutual separation at all depths, one in which the separation increases with depth, but their size 
is still unchanged, and one in which both size and separation increase with depth, which is the standard flow 
structure. It is concluded that the third possibility with fewer and thicker downdrafts in deeper layers remains 
the most plausible, but it may be unable to explain the suspected absence of large-scale flows with speeds and 
scales expected from MLT. 

Subject headings: convection — turbulence — Sun: granulation 


1. INTRODUCTION 

Late-type stars such as our Sun have outer convection 
zones. The observed solar granulation is a surface manifes¬ 
tation of their existence. Solar granulation and the first few 
megameters (Mm) of the convection zone have been modeled 
successfully using mixing length theory (MLT) and numerical 
simulations with realistic physics included (Stein & Nordlund 
1989, 1998; Vogler et al. 2005; Gudiksen et al. 2011; Ereytag 
et al. 2012). As a function of depth, the simulations repro¬ 
duce some essential features predicted by MLT, in particular 
the depth dependence of the turbulent root mean square (rms) 
velocity, Urms ~ {Fenth/py^^, where i^enth is the enthalpy 
flux and p is the density. Simulations also seem to confirm an 
important hypothesis of MLT regarding the gradual increase 
of the typical convective time and length scales with depth. 
Given such agreements, there was never any reason to ques¬ 
tion our basic understanding of convection. 

In recent years, local helioseismology has allowed us to de¬ 
termine subsurface flow velocities (Duvall et al. 1997; Gizon 
& Birch 2005). Helioseismic observations by Hanasoge et 
al. (2010, 2012) using the deep-focusing time-distance tech¬ 
nique have not, however, detected large-scale convection ve¬ 
locities at the expected levels (Miesch et al. 2012); see also 
Gizon & Birch (2012) for a comparison with both global sim¬ 
ulations and radiation hydrodynamics simulations with real¬ 
istic radiation and ionization physics included. Greer et al. 
(2015) have suggested that the approach of Hanasoge et al. 
(2012) may remove too much signal over the time span of the 
measurements. Using instead ring diagram analysis (Gizon 
&. Birch 2005) with appropriately assembled averaging ker¬ 
nels to focus on the deeper layers, Greer et al. (2015) instead 
And values of the turbulent rms velocities that are consistent 
with conventional wisdom. Moreover, they And that at large 


length scales corresponding to spherical harmonic degrees of 
30 or below (corresponding to scales of about 140 Mm and 
larger), the rms velocity actually increases with depth. This 
in itself is remarkable, because velocity perturbations from 
deeper layers were expected to be transmitted to the surface 
almost unimpededly (Stix 1981), unlike perturbations in the 
heat flux, which are efficiently being screened by the convec¬ 
tion (Spruit 1977). Thus, an increase in horizontal velocity 
power with depth is unexpected. However, van Ballegooijen 
(1986) pointed out that for the scale and depth of giant cells, 
the screening might be large enough to allow giant cell con¬ 
vection of 100ms“^ to result in only 10ms“^ at the surface, 
which would be compatible with observations (Hathaway et 
al. 2013). The presence of the near-surface shear layer of the 
Sun (Schou et al. 1998), might enhance the screening further. 

Given that the deeply focused kernels used by Greer et al. 
(2015) can still have 5-10% sensitivity to Doppler shifts aris¬ 
ing from flows in the upper layers, such near-surface flows 
could still leave an imprint on the signal. This effect would 
be exaggerated further if the near-surface flows were stronger 
compared to those in deeper layers, as was theoretically ex¬ 
pected. On the other hand, if convection in the surface layers 
is of smaller scale, the signal from those layers will to a large 
extent be averaged out despite its larger amplitude. It would 
obviously be important to examine this more thoroughly by 
applying the kernels of Greer et al. (2015) to realistic simu¬ 
lations. Unfortunately, this has not yet been attempted. Con¬ 
versely, the deep-focusing technique of Hanasoge et al. (2012) 
could be extended to allow for imaging of the deeper flow 
structures and thereby a direct comparison of individual tur¬ 
bulent eddies with those detected by Greer et al. (2015). 

It should be mentioned that small flow speeds of giant cell 
convection have been found by correlation tracking of super- 
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granule proper motions (Hathaway et al. 2013). The typical 
velocities are of the order of 20ms“^ at spherical degrees 
of around 10. Those speeds are still an order of magnitude 
above the helioseismic upper limits of Hanasoge et al. (2012), 
but Hanasoge & Sreenivasan (2014) argue that these surface 
measurements translate to lower flow speeds in the deeper and 
denser layers when considering mass conservation. This ar¬ 
gument may be too naive and would obviously be in conflict 
with the results of Greer et al. (2015), which show an increase 
with depth at low spherical harmonic degrees. On the other 
hand, Bogart et al. (2015) have argued that the flows reported 
by Hathaway et al. (2013) may be non-convective in nature 
and, in fact, magnetically driven and perhaps related to the 
torsional oscillations. 

Quantitatively, the horizontal flow speed at different scales 
is characterized by the kinetic energy spectrum, E{k), where 
k is the wavenumber or inverse length scale. Stein et al. 
(2007) have shown that the spectra of surface Dopplergrams 
and correlation tracking collapse with those of Stein & Nord- 
lund onto a single graph such that the spectral velocity 
[kE{k)Y/^ is proportional to k, i.e., E{k) oc k; see also 
Nordlund et al. (2009). This is a remarkable agreement be¬ 
tween simulations and observations. While the origin of such 
a spectrum is theoretically not understood, it should be em¬ 
phasized that these statements concern the horizontal veloc¬ 
ity at the surface and do not address the controversy regard¬ 
ing the spectrum at larger depths, where the flow speeds at 
large length scales may still be either larger (e.g. Greer et al. 
2015) or smaller (e.g. Featherstone & Hindman 2016) than 
those at the surface. If most of the kinetic energy were to 
reside on small scales throughout the convection zone even 
at a depth of several tens of megameters and beneath, E{k) 
is expected to decrease toward smaller wavenumbers k either 
like white noise oc k^ or maybe even with a Batchelor spec¬ 
trum oc (Davidson 2004). However, even in simulations 
of forced isothermally stratified turbulence, in which there is 
no larger-scale driving from thermal buoyancy, there is a shal¬ 
lower spectrum, E{k) oc (Losada et al. 2013). By con¬ 
trast, if one imagines the flow to be anelastic so that the mass 
flux pit = V X i/) can be written as a vectorial stream func¬ 
tion i/), and if t/; is given by white noise, one should expect a 
Batchelor spectrum. This is qualitatively similar to the re¬ 
sults of granule tracking, which reveal an intermediate scaling 
proportional to k^ (Rieutord et al. 2008; Roudier et al. 2012). 
While the steeper k^ spectrum for k < k{ might leave some 
hope that the results of Hanasoge et al. (2012) could be recon¬ 
ciled with these theoretical and observational constraints, this 
would be virtually impossible for the linear energy spectrum, 
E{k) oc k. 

The issue has recently been examined by Lord et al. (2014), 
who have shown that simulation results with the MURaM 
code (Vogler et al. 2005) can be reproduced with a model 
that is composed of a continuous hierarchy of layers, each 
with its own driving scale at a scale of four local density scale 
heights. At wavenumbers below that driving scale, the spec¬ 
tral power falls off in a certain way. The horizontal surface 
spectrum is a superposition of these contributions, each of 
which is assumed to decay with height, contributing therefore 
progressively less with depth. However, both their model and 
the MURaM results (Lord 2014) show an order of magnitude 
more power (factor 2.5 in rms velocity) than the Sun at small 
wavenumbers. To understand this, they show that the obser¬ 
vations can be reproduced if below a depth of 10 Mm the con¬ 


vective energy flux is augmented by an artificially added flux 
term, which thus significantly reduces the rms velocity of the 
resolved flow field. They speculate that the flow-temperature 
correlation entering the enthalpy flux may be larger in the Sun, 
possibly being caused by a magnetic held that may maintain 
flow correlations and boost the transport of convective flux by 
smaller scales. This has been partially confirmed by Hotta et 
al. (2015), who have discussed the possible role of small-scale 
magnetic fields in suppressing the formation of large-scale 
flows. Furthermore, Featherstone & Hindman (2016) have 
found that with increasing Rayleigh numbers, there is more 
kinetic energy at small scales and less at large scales such that 
the total kinetic energy is unchanged by this rearrangement of 
energy. 

Given the importance of this subject, it is worthwhile re¬ 
viewing possible shortcomings in our theoretical understand¬ 
ing and numerical modeling of stellar convection. Both global 
and local surface simulations of solar-type convection would 
become numerically unstable with just the physical values of 
viscosity and radiative diffusivity. Therefore, the radiative 
flux (in the optically thick layers) is modified in one of two 
possible ways, (i) The contribution from temperature fluctua¬ 
tions is greatly enhanced, i.e., 

F,,d = -KVT ^ -iTVT-iTsGsV(T-T), (1) 

where T and T are the actual and horizontally averaged tem¬ 
peratures, respectively, K is the radiative conductivity, and 
Ksgs is a subgrid scale (SGS) conductivity. The latter is 
enhanced by many orders of magnitude relative to K. Fur¬ 
thermore, numerical diffusion operators often do not translate 
in any obvious way to the physical operators, (ii) Alterna¬ 
tively, in direct numerical simulations (DNS), one uses physi¬ 
cal viscosity and diffusivity operators, i.e., the replacement in 
Equation (1) is not invoked, but the coefficients are enhanced 
{K iTenh) and exceed the physical ones by many orders of 
magnitude. Both approaches are problematic. 

In many global DNS with enhanced coefficients (e.g., 
Kapyla et al. 2013), the lower boundary is closed, so at the 
bottom of the domain all the energy is carried by radiation 
alone. By choosing an enhanced radiative diffusivity, the 
radiative flux is increased by a corresponding amount, and 
therefore also the total flux. The luminosity in those simula¬ 
tions can exceed the solar value by several orders of magni¬ 
tude. This has a series of consequences. Most notably, the 
convective velocities are too high in the upper parts where 
most of the flux is carried by convection (Kapyla et al. 2013). 
There are two ways to avoid this problem. One is to use 
simulations with a polytropic hydrostatic reference solution, 
whose polytropic index is close enough to the adiabatic one 
so that the convective flux is everywhere a small fraction of 
the radiative one (Brandenburg et al. 2005), which reduces 
the convective velocities. Alternatively, one applies the en¬ 
hanced radiative conductivity only to the temperature fluctua¬ 
tions so as not to disturb the very small radiative energy flux, 
—A'VT, compared with — iTenhVT, which it would have 
been in the DNS approach without subtracting T. The tem¬ 
perature smoothing implied by invoking alternative (i) is also 
necessary in the simulations with realistic opacities, because 
the Peclet number (which is similar to the Reynolds number) 
based on the physical value of K would reach values above 
10^°, which cannot be handled by a DNS (Barekat &. Bran¬ 
denburg 2014). In the global simulations, it is common to 
use the specific entropy gradient instead of V(T — T) (e.g. 
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Kapyla et al. 2013). This is equivalent if T is close to the 
adiabatic value. However, the diffusion coefficient in the SGS 
term can easily be hve orders of magnitude larger than the 
physical one acting on the mean stratihcation. This has the 
consequence of suppressing small-scale turbulent flows and 
entropy fluctuations, which may prematurely damp the low- 
entropy fluid parcels that originate within the strong cooling 
layer at the surface and which will be discussed as entropy 
rain in the bulk of this paper. If such low-entropy contri¬ 
butions are poorly captured by the simulation, this could be 
compensated for by a sufficiently strong contribution propor¬ 
tional to the superadiabatic gradient, which in turn would give 
rise to the excitation of flows on much larger scales than what 
would be compatible with the observations discussed by Lord 
et al. (2014). The sensitivity of large-scale flow excitation 
to small changes in the superadiabatic gradient was also dis¬ 
cussed by Cossette & Rast (2016), who studied models with 
different types of surface driving. It is further supported by 
the work of Featherstone & Hindman (2016), who found that 
there is more kinetic energy at small scales and less at large 
scales as the Rayleigh number increases. 

There is yet another problem that concerns the global simu¬ 
lations in which a predetermined profile K{z) is used instead 
of calculating its local value with a physical opacity. Such 
an approach has been used routinely in studies of compress¬ 
ible convection, especially when stably and unstably stratified 
layers are combined (Hurlburt et al. 1986; Brandenburg et al. 

1996) . The simulations of Kapyla et al. (2013) adopted a pro¬ 
file that yields a negative (unstable) radial entropy gradient 
through most of the layer (see the inset of their Fig. 3). How¬ 
ever, as will be pointed out below, it is only a tiny surface 
layer in which the non-convective model is unstable. The 
rest of the model is a priori stable to convection, but it be¬ 
comes marginally stable (or perhaps slightly unstable) as a 
consequence of the resulting turbulence leading to bulk mix¬ 
ing across the deeper layers. This is quite contrary to the mod¬ 
els with a predefined K{z) profile, which would be unstable 
by construction over the entire depth of the convection zone. 

We should emphasize that the non-convective solution is 
mainly of academic interest. It is used to compute, for exam¬ 
ple, the Rayleigh number, a measure of the degree of insta¬ 
bility. However, there is no doubt that the actual stratification 
of the Sun is close to marginal stability down to a fractional 
radius of 0.71 before turning decidedly stable, as confirmed 
by helioseismology (Christensen-Dalsgaard et al. 1991; Basu 

1997) . Nevertheless, the concept of convection being driven 
by surface cooling rather than heating from below has been 
promoted in a number of papers by Stein & Nordlund (1989, 

1998) and elaborated upon by Spruit (1997), who introduced 
the idea that flows in the deeper parts are being driven nonlo- 
cally through the entropy rain from the surface. The question 
is to what extent this affects our understanding of the speed 
and especially the typical scales of convective motions, how 
MLT models would need to be modified, and whether this 
might have any bearing on the interpretation of the observed 
flow amplitudes at large scales. 

To include the nonlocal effects described by Spruit (1997), 
we must look for a contribution to the enthalpy flux that is 
not related to the local entropy gradient.' Such a term has 

* As will be discussed in Section 6, Spruit (1997) argues that the low- 
entropy material from the top always leads to a negative mean entropy gra¬ 
dient. However, here we argue that the mean entropy gradient is only one 
contribution to a mean-field (here one-dimensional) parameterization of the 


been identified in the meteorological context by Deardorff 
(1966), who describes it as a counter-gradient flux. In Dear¬ 
dorff (1972), he derives an expression for this flux, which de¬ 
pends on the local temperature or entropy fluctuations which, 
in his case, come from measurements. In the present case, we 
assume that such fluctuations have their origin in what Spruit 
(1997) refers to as threads, which are thin downdrafts on an 
almost perfectly isentropic passive fluid upflow between the 
threads. Spruit already emphasized back then (p. 406) the dy¬ 
namical importance of the “forest of narrow cool threads that 
are produced at the solar surface.” He further suggested (p. 
406) that their absence in simulations with a fixed top bound¬ 
ary at some depth below the actual surface might be “the rea¬ 
son why they would produce large-scale flows with ampli¬ 
tudes that are about two orders of magnitude larger than ob¬ 
served.” 

In the scenario described by Spruit (1997), convection is 
driven solely within the surface layers. This is superficially 
reminiscent of convective overshoot, as modeled in some of 
the aforementioned papers (Hurlburt et al. 1986; Brandenburg 
et al. 1996), where convective flows are driven into stable lay¬ 
ers. The extent of overshoot depends not only on the stiffness 
of the stable gradient that the convective plumes are flowing 
into (Hurlburt et al. 1994), but also on the flow speed. The 
depth of such a layer can therefore not be determined a priori 
from purely hydrostatic stability considerations. The lower 
boundary might then not be sharp. As explained above, in a 
hydrostatic non-convecting reference model, only a very thin 
surface layer is a priori unstable to convection. The rest is 
made unstable purely by bulk mixing. If entropy rain con¬ 
vection is a nonlocal phenomenon, the extent of convection 
should depend on surface properties and cannot be predicted 
from the local entropy gradient, similarly to convective over¬ 
shoot. One might then not be able to understand the relatively 
sharp demarcation at the bottom of the convection zone, as 
found in global helioseismology. Of course, once convec¬ 
tion has become fully developed, the low entropy elements 
descend into buoyantly neutral layers, which is quite different 
from the usual overshoot. Furthermore, the usual overshoot 
layer is characterized by negative buoyancy and therefore a 
downward enthalpy flux (Hurlburt et al. 1986). An impor¬ 
tant purpose of the present paper is to produce a quantitative 
model and to demonstrate that, with a hypothetical nonlocal 
contribution to the flux, it is possible to obtain models that 
still have a sharply defined lower boundary. 

A quantitative model of convection with entropy rain, even 
with the somewhat hypothetical Deardorff term included, 
would also be useful to illustrate the qualitative nature of the 
resulting stratification, and to show whether it is weakly super 
or subadiabatic. Indeed, as already argued by Spruit (1997), if 
the local mass fraction of entropy deficient material, descend¬ 
ing from the cooling surface, decreases with depth and if the 
stratification outside the entropy rain were exactly isentropic, 
the resulting horizontally averaged entropy would increase 
with depth. This suggests that the entropy rain itself can make 
an otherwise vanishing radial entropy gradient negative and 
therefore Schwarzschild unstable, as seen in surface simu¬ 
lations. This raises two important questions. First, to what 
extent is such a stratification affected by radiative heating— 
especially toward the bottom of the convection zone, where it 
would tend to produce a positive (stable) mean entropy gradi- 

enthalpy flux, and that there is another one that is not locally connected with 
it. 
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ent in the upflows. Second, would such a negative (unstable) 
mean entropy gradient always lead to giant cell convection, as 
has been seen in global convection simulations (Miesch et al. 
2008), or could the stratification still be stable to the excitation 
of large-scale flows? 

We postpone addressing the two aforementioned issues 
connected with the qualitative reasoning of Spruit (1997) to 
the end of the paper (Section 6), and begin by highlighting the 
less commonly known fact that in the non-convecting refer¬ 
ence model, only the layer in the top 1 Mm is convectively un¬ 
stable (Section 2). We then explain the nature of the Deardorff 
flux (Section 3), present a correspondingly modified mixing 
length model (Section 4) and give some illustrative numerical 
solutions (Section 5). We discuss alternative explanations for 
the lack of giant cell convection in Section 6 and present our 
conclusions in Section 7. 

2. A HIGHLY UNSTABLE SURFACE LAYER 

The main argument for the existence of a highly unstable 
surface layer comes from the consideration of the associated 
non-convective reference solution, where the flux is forced 
to be transported by radiation only, i.e. F = Frad- This is 
something that is not normally considered in stellar physics, 
because we know that such a solution would be unstable to 
convection and would therefore never be realized. However, 
as mentioned in the introduction, this solution is of certain 
academic interest. We postpone the discussion of an explicit 
numerical solution to Section 5 and present in this section the 
basic argument only at a qualitative level, making reference 
to earlier numerical calculations by Barekat & Brandenburg 
(2014) for a simple model. They considered an opacity law of 
the form 

K = Ato(p/po)“(T/To)^ ( 2 ) 

where a and b are adjustable parameters, po and Tq are ref¬ 
erence values for density and temperature, respectively, and 
kq gives the overall magnitude of the opacity. The essential 
point to note here is that the exponents a and b determine the 
gradient of specific entropy in a purely non-convecting refer¬ 
ence model. In thermodynamic equilibrium, the radiative flux 
must be constant, i.e., 

Frad = —K dT/dz = const, (3) 

where K = 16crsBF^/(3Kp) is the radiative conductivity 
with (TsB being the Stefan-Boltzmann constant, and z is the 
vertical coordinate in a Cartesian coordinate system. For the 
simple opacity law ( 2 ), but with 

6 < 4 -b a, (4) 

the optically thick regime is characterized by a constant tem¬ 
perature gradient and therefore K = const (see Appendix A). 
We have then a polytropic stratification with p oc T", where 

n = (3 - &)/(! -b a) (5) 

is the polytropic index. It is larger than —1 when Equation (4) 
is obeyed, i.e., when the pressure decreases with height, as 
is required for a physically meaningful solution. For a ra¬ 
tio of specific heats of 7 = 5/3, the value of n for marginal 
Schwarzschild stability is ricrit = 1 /( 7 ~ 1 ) = 3/2. Inmost of 
the solar convection zone the dominant opacity is the bound- 
free absorption owing to the absorption of light during ion¬ 
ization of a bound electron, which is well described by the 
Kramers-type opacity law with a = 1 and b = —11% so 


n = 3.25, which corresponds to a Schwarzschild-stable solu¬ 
tion. Only near the surface, at temperatures typically below 
15,000 K, the dominant opacity is the H“ opacity. It can no 
longer be approximated by a simple power law of the type 
given by Equation (2). However, in limited density and tem¬ 
perature ranges, certain values of a and b can tentatively be 
specified, e.g., a = 0.5 and b = 7...18. Clearly, for all these 
values the constraint (4) is violated, so the hydrostatic strati¬ 
fication is no longer polytropic, but solutions can still be con¬ 
structed numerically (Barekat & Brandenburg 2014) and they 
demonstrate, not surprisingly, that the stratification is highly 
unstable. What is more surprising is the fact that even with a 
combined opacity law of the form 

K = “b ; (6) 

where KKr and kh- are given by Equation (2) with suitable 
exponents a and b, the solutions of the non-convective ref¬ 
erence state is unstable only over a depth of approximately 
1 Mm. We return to this at the end of Section 5, where we 
present numerical solutions. 

Of course, as stated earlier, the non-convective reference 
state is only of academic interest. Already with standard 
MET (Biermann 1932; Vitense 1953), which allows for a non¬ 
vanishing enthalpy flux, one finds a vastly extended convec¬ 
tion zone with a depth of the order of 100 Mm (Biermann 
1938). However, the question now is how this can be affected 
by the presence of the Deardorff flux. This will be the subject 
of the rest of this paper. If the Deardorff flux were to become 
dominant and the stratification subadiabatic, it would become 
locally stable to the onset of convection. One might further 
speculate that the typical scale would no longer be controlled 
by the local pressure scale height, but it might be imprinted 
from the downdraft pattern just beneath the surface and be 
therefore comparable to the granulation scale or at least the 
supergranulation scale (Cossette & Rast 2016). This can have 
other important consequences that will also be addressed in 
this paper. It should be pointed out, however, that the sur¬ 
face simulations have so far not produced evidence for sub¬ 
adiabatic stratification. On the other hand, the models pre¬ 
sented below predict subadiabatic stratification only a certain 
distance below the surface, depending on ill-known input pa¬ 
rameters. 

3. THE DEARDORFF FLUX 
3.1. Derivation 

In the meteorological context, counter-gradient heat flux 
terms have been noticed for a long time (Ertel 1942; Priestley 
& Swinbank 1947; Deardorff 1966). They appear naturally 
when calculating the enthalpy flux Fenth using the r approxi¬ 
mation in its minimalistic form (e.g. Blackman & Field 2003). 
In this approach, one computes the time derivative of Fenth- 
In the absence of ionization effects, Fenth can be written as 

Fenth = pU^CpT, (7) 

where the overbar denotes horizontal averaging. In standard 
MET, the enthalpy flux is usually referred to as the convec¬ 
tive flux and the kinetic energy flux vanishes because of the 
assumed perfect symmetry between up- and downflows. In 
deeper layers especially, however, this is not justified (Stein 
et al. 2009), so this restriction will later be relaxed. 

In the case of strongly stratified layers, it is convenient to 
use pressure P and specific entropy S as thermodynamic vari¬ 
ables, because we later want to ignore pressure fluctuations 
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on the grounds that pressure disturbances are quickly equi¬ 
librated by sound waves. Furthermore, we will restrict our¬ 
selves to second order correlations in Equation (7) and thus 
to fluctuations only in the correlation of specific entropy and 
velocity. Up to some reference value, we have 

5/cp = In T-Vad In P, (8) 

where Vad = 1 — 1 / 7 , and 7 = cp/cy is the ratio of specific 
heats at constant pressure and constant volume, respectively, 
all constants for a perfect gas with a fixed degree of ionization. 
Ignoring pressure variations and linearizing i5 In T = ST/T, 
we can replace cpST by TSS. In the following, we denote 
fluctuations by lower case characters, i.e., S = S + s, where 
s = 5S are used interchangeably. As argued above, other fluc¬ 
tuations are omitted. We also ignore mean flows (U = 0), so 
there are only velocity fluctuations {U = u). Focusing thus 
on the dominant contribution proportional to ujs, we have 

Tenth =pTvTs. (9) 

Next, we write the time derivative of Fenth as 

dFenth/dt = pT (u^s + iizs) , (10) 

where dots denote partial time derivatives and changes of the 
background state have been neglected. Using the governing 
equations for Ui and s, we have (see Appendix B) 

S = -UjVjS - s/TcooI---, ( 11 ) 

Ui = -gisjcp + (12) 

where gi is the zth component of the gravitational acceleration 
in Cartesian coordinates, namely g = (0,0, —g), the ellipses 
refer to terms that are nonlinear in the fluctuations. 


(15) have been verified numerically for passive scalar trans¬ 
port (Brandenburg et al. 2004). Among other things, they 
found that the time derivative on the left-hand side of Equa¬ 
tion (14) restores causality by turning the otherwise parabolic 
heat equation into a hyperbolic wave equation. Here, how¬ 
ever, we are interested in slow variations such that the time 
derivative of Fenth can be neglected. (Note, however, that this 
does not imply that the cooling term will be neglected.) 

Since the relaxation term in Equation (15) is similar to the 
cooling term in Equation (14), we can combine the two by 
introducing the reduced relaxation time Tred defined through 


We can then solve for Fgnthj which appears on the right-hand 
sides (rhs) of Equations (14) and (15), and find Fenth = Fg+ 
Fd, where 

FG = -^T^edu'^msP'r'^S, (17) 

FB=-Treds‘^ gpT/cp (18) 

are the ordinary gradient and the new Deardorff fluxes, re¬ 
spectively, and anisotropies have been ignored for the benefit 
of simpler notation. Thus, we write UiUj « ^SijU^^^, where 
Urms is the rms velocity of the turbulence. Equations (17) 
and (18) are written in vectorial forms to highlight the direc¬ 
tions of the fluxes; Fq is counter-gradient and Fd is coun¬ 
tergravity. In Equation (17), the term IrredU^ms = Xt is the 
turbulent thermal diffusivity. 

In view of astrophysical applications, we replace the spe¬ 
cific entropy gradient by the commonly defined superadia- 
batic gradient, i.e., 

-d(S/cp)/dz = (V - Vad)/Fp, (19) 


Teooi = I'C-ykf with L = ikf/{3 + kf) (13) 

is the inverse heating and cooling time owing to radiation 

(Unno & Spiegel 1966; Edwards 1990), Cj = WctsbT^/ pcp 
is the photon diffusion speed (Barekat & Brandenburg 2014), 
i = 1 /kp is the photon mean-free path, and k{ is the typi¬ 
cal wavenumber of the fluctuations, which may be associated 
with the inverse mixing length used routinely in MET. We 
mention at this point that the nonlocal nature of the entropy 
rain must lead to yet another contribution in Equation (11^ 
We expect this to be the result of the nonlinear term (indi¬ 
cated by the ellipses), of which a part later gives rise to the 
Vd term. This will be motivated further at the end of this sec¬ 
tion and in Section 3.3, where it will be included in the final 
expression for Vd- 

Using Equations (11) and (12) in Equation (10), we have 
— h—=pT (-upIj'^jS - gis'^/cp) - - - \-Ti, (14) 

ut \ / 'T’cool 

where the term is the Deardorff flux and the 7) refer to 
triple correlations that will be approximated by the quadratic 
correlation Fenth as 

T, = -Fr^^lr, (15) 

where r is a relaxation time due to turbulence (Blackman & 
Field 2003), which will later be identified with the turnover 
time. This procedure, in which correlations with pressure 
fluctuations are also neglected, is called the minimal r ap¬ 
proximation. Important aspects of the closure assumption 


where V = d In T/d In F is the double-logarithmic tempera¬ 
ture gradient and 7/p = —(dlnF/dz)”^ is the pressure scale 
height. Thus, we arrive at 

Fenth = |pCpT {TredU^ms/Hp) ~ ^ad + Vd) , (20) 

where Vd is a new contribution to standard MET, which re¬ 
sults from Fd. Using the terms on the rhs of Equation (18), 
we can write it explicitly as 

Vd = (3/7) (^/ 4 ) Ma- 2 , (21) 

where Ma = Urms/cs is the Mach number of the turbulence 
and Cg is the sound speed with Cg = 7p7/p. Note that this 
Mach number dependence arises in order to cancel the cor¬ 
responding u^jns factor in the definition (20). The Deardorff 
term is a contribution to the flux that is always directed out¬ 
ward and results from the transport of fluid elements with en¬ 
tropy fluctuations of either sign, as is illustrated in Figure 1. 

3.2. Physical interpretation 

We recall that in a stably (unstably) stratified layer, shown 
in Figure 1(a) and (b), and under the assumption of pressure 
equilibrium, the entropy perturbation of a blob with respect 
to its current surroundings decreases (increases) after a small 
ascent. By contrast, in the marginally or nearly marginally 
stratified cases, the perturbation remains constant if one ig¬ 
nores mixing with the surroundings; see panels (c) and (d). 
Therefore, a positive entropy perturbation (s > 0) always 
corresponds to a positive temperature perturbation and a neg¬ 
ative density perturbation. Thus, the fluid parcel is buoyant 
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Fig. 1.— Sketch illustrating overshoot in a stably stratified layer (a), growth 
of perturbations in an unstably stratified layer (b), buoyant rise of a blob with 
positive entropy perturbation (c), and descent of a blob with negative entropy 
perturbation and hence negative buoyancy (d). In cases (b)~(d), the turbulent 
flux is upward (increasing z). Case (d) is relevant to entropy rain. 


and moves upward {Uz > 0), so UzS > 0, giving a positive 
contribution to i^enth- Likewise, for s < 0, the parcel is cooler 
and heavier and sinks (uz < 0), and again, UzS > 0, giving a 
positive contribution to -Fenth- 

In laboratory convection, strong positive entropy perturba¬ 
tions can be driven at the lower boundary, and strong nega¬ 
tive ones at the top boundary. In the Sun, however, only the 
strong radiative cooling in the photosphere provides a signifi¬ 
cant source of entropy perturbations. This became clear with 
the emergence of the realistic solar convection simulations of 
Stein & Nordlund (1989). These simulations showed for the 
first time the strong vertical asymmetry of solar convection 
resulting from the physics of ionization and strong cooling 
via the H” opacity. This led to the notion of entropy rain 
and Spruit’s description of solar convection as a nonlocal phe¬ 
nomenon that is driven solely by surface cooling. Our consid¬ 
erations suggest that the resulting (one-dimensional) mean- 
field energy flux is parameterized not just in terms of the local 
superadiabatic gradient. 

In analogy with laboratory convection (Heslot et al. 1987; 
Castaing et al. 1989), Spruit refers to the downdrafts as 
threads. Interestingly, the laboratory experiments show that 
these threads persist even at large Rayleigh numbers. For such 
threads to persist, one could imagine them to be stabilized 
by their intrinsic vorticity, akin to a Hill vortex (Hill 1894). 
These are vortex rings that can stay concentrated over large 
distances. Vortex rings have also been reported by Stein et 
al. (2009), but many of those are associated with mushroom 
shapes, indicating that they widen and soon break up and mix 
with their surroundings. Numerical studies (see Appendix C) 
suggest that Hill vortices also survive in a highly stratified 
isothermal layer, and that their persistence increases signifi¬ 
cantly with resolution and decreasing viscosity. On the other 


hand, the presence of background turbulence provides an ef¬ 
fective turbulent viscosity, which contributes to a premature 
breakup of small-scale vortices. However, those studies were 
done with an isothermal equation of state and thus ignore the 
effects of continued driving by a persistent negative entropy 
perturbation between the vortex and its surroundings. 

Studies of the Deardorff flux in the meteorological context 
suggest that the flux-carrying plumes are associated with so- 
called coherent structures (De Roode et al. 2004) and those 
may not simply be the vortex-like structures envisaged above. 
It is clear, however, that the structures seen in the Earth’s at¬ 
mosphere are driven by the boundary layer on the ground, 
and sometimes also by the upper inversion layer (De Roode 
et al. 2004). Pleim (2007) discusses the Deardorff flux in con¬ 
nection with other parameterizations of nonlocal fluxes such 
as the transilient matrix approach (Stull 1984, 1993). The 
close connection between counter-gradient fluxes and nonlo¬ 
cal transport has been elaborated upon by van Dop &. Verver 
(2001) and Buske et al. (2007). Thus, while the approach pre¬ 
sented in Section 3.1 may capture the physical phenomenon 
of the Deardorff flux qualitatively correctly, it is quite possi¬ 
ble that the nature of the underlying coherent structures may 
require additional refinements. 

3.3. Depth dependence of the Deardorff term 

To estimate the resulting depth dependence of the Deardorff 
term in Equation (21), we must know how the associated 
with the Deardorff term varies with depth. If the entropy rain 
was purely of the form of Hill vortex-like structures, as dis¬ 
cussed above, their filling factor fg would decrease with in¬ 
creasing depth and density like 

/. ( 22 ) 

where C, = 0.8 has been found for spherical vortex struc¬ 
tures descending in an isothermally stratified layer; see Ap¬ 
pendix C. However, Equation (23) is not contingent on Hill 
vortices, but is a consequence of downdrafts along a density 
gradient. Eor purely spherical compression one would expect 
C, = 2/3, while for horizontal compression one has C = 1- 

If we neglect non-ideal (radiative or viscous) effects, as 
well as entrainment between up- and downflows, the differ¬ 
ence AF = Sy — Si between the entropy of the slowly ris¬ 
ing surroundings, Sy, and that in the downward propagating 
vortex. Si, would remain constant and equal to the entropy 
deficit AFo suffered by overturning motions at the surface, 
which is the only location where radiative losses are signif¬ 
icant. On sufficiently short length scales, however, radiative 
heating from the surroundings would erode this entropy dif¬ 
ference and lead to a decrease of the effective /S.S. We model 
this by considering C, an adjustable parameter that is increased 
relative to the value 0.8 that we expect in the absence of ra¬ 
diation, i.e., fs decreases more strongly. The resulting frac¬ 
tional entropy difference, fg ASq, therefore decreases faster 
with depth than for ( = 0.8. 

As will become clear from the considerations below, if there 
is no entrainment from the upflows into the downdrafts and all 
of the downflows were confined to the small surface area of 
the vortex structures, the resulting downward-directed kinetic 
energy flux would increase with depth and eventually exceed 
the enthalpy flux. This would be unphysical, because con¬ 
vection should lead to outward energy transport. However, if 
there is entrainment with associated mixing, the actual filling 
factor (fractional area) of downflows, which we denote by /, 
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would be larger and the mean entropy difference /S.S would 
be diluted correspondingly such that 

fAS = fsASo (ocp-^AS-o). (23) 

We recall that non-ideal effects can be captured by choosing 
C, > 0.8 in the prescription (23). The physics of entrainment 
has been modeled in the stellar context by Rieutord & Zahn 
(1995) and the extent of entrainment has been quantified in the 
realistic surface simulations by Trampedach & Stein (2011), 
who determined the entrainment length scale as the typical 
scale height of the mass flux in the up- and downflows sepa¬ 
rately. They call this scale the mass mixing length and And its 
value to be comparable to Hp. 

To estimate the depth dependence of various horizontal av¬ 
erages, and to compute enthalpy and kinetic energy fluxes, 
we now compute various averages in the two-stream approx¬ 
imation, in which all horizontal averages are the sum of the 
fraction 1 — / of the value in upflows and the fraction / of 
the value in downflows. Thus, for the mean specific entropy 
stratification we have 

S={1- f)S^ + fS^ = S^- fAS, (24) 

where and Si are the mean specific entropies in up- and 
downflows, and AS = S^ — Si is their difference. We expect 
that Si ~ Si will be constant if there is no heating in the 

upwellings, while Si fa {1 - fs/f)Si + {fs/DSi will be 
dominated by contributions from Si due to entrainment. 

To calculate s^, we must first compute the fluctuating quan¬ 
tity s = S — S and then average the squared values in up- and 
downflows, which, using Equation (24), gives 

^ = (1 - /) (^t -Sf + f (Si - sf = f {ASf, (25) 

where / = (1 — /)/ has been introduced as a shorthand. 
Analogously to the specific entropy, we And for the velocity 

C7, = (1 - f)Ui + fUi = Ui-f AU, (26) 

where Ui > 0 and Ui < 0 are the mean up- and downflow 
velocities with A[/ = Ui — Ui- Since the densities in up- and 
downflows are nearly the same, especially in deeper layers, 
we have Uz ~ Q from mass conservation. With this, we And 

= ^z = (^-f)u\ + ful = f (A[7)^ (27) 

and thus C/t/urms = [//(1 — /)]^^^- Similarly, we calculate 

^ = (1 -/)[/' + ful = -/(I - 2/) (AC7)3, (28) 

which is negative for / < 1/2, and Anally 

uls = (1 - /) C7t5t + / UiSi = f AU AS. (29) 

Thus, the kinetic energy flux, pu'^Uz/2, which reduces to 
p u1/2 in this two-stream approximation, is given by 

Tkin — 0kin Perms’ (30) 

where i/kin = (1/2 — /)//^^^ is a positive prefactor (corre¬ 
sponding to downward kinetic energy flux) if / < 1/2. Stein 
et al. (2009) And / Ri 1/3, nearly independently of depth, 
which yields i/kin ~ \/2/4 r; 0.35; see Table 1, where we 
list 0kin and -Ui/urms = [(1 - for selected val¬ 

ues of /. The enthalpy flux is proportional to and, using 


_ _ TABLE 1 

4‘kin, —Ui/Uruis, AND 17.|-/ltrms FOR SELECTED VALUES OE /. 

/ _ 1/2 1/3 0.14 0.015 0.0006 

_((>kin 0 0.35 1 4 20 

^4,^rms 1 1.4 2.5 8 40 

Ui/Uru,s 1 0.7 0.4 0.12 0.025 

Equations (9) and (29) together with Equations (25) and (27), 
we And E’enth = 7>T UrmsSrms- In Appendix D we show that 
Tsrms = u?ms^ff^p/(aMLTVad), where omlt = 1/8 is a 
geometric factor in standard MLT^ This leads to 

Eenth — /^enth PTlj-jns (31) 

with (/enth = fcfEp/(aMLTVad)- TWs yields /ignth ~ 20, 
which is rather large. By contrast, Brandenburg et al. (2005) 

determined a quantity such that (penth = fa 4. 

We see that the presence of a kinetic energy flux just mod¬ 
ifies the usual expression for the convective flux, which then 
becomes the sum of enthalpy and kinetic energy fluxes; see 
Section 3.1. This is compatible with recent simulations of 
R. E. Stein (2016, private communication), in which the frac¬ 
tional kinetic energy flux increaes toward the deeper parts 
(^ 40 Mm) of the domain. 

When / becomes small (< 0.14), /)kin exceeds unity and 
for / < 0.015, (/kin exceeds the estimate cpenth ~ 4 found by 
Brandenburg et al. (2005), so the sum of enthalpy and kinetic 
energy fluxes may become negative, which appears unphysi¬ 
cal. In that case, the idea of reconciling the results of Hana- 
soge et al. (2012) with such convection transporting the solar 
luminosity could be problematic, unless both up- and down¬ 
flows were to occur on sufficiently small scales. In the case 
/ = 0.14, the resulting Ui would still not be particularly fast 
and only 2.5 times larger than Urms; for / = 0.015 we have 
— Ui fa 8urms; see Table 1. Conversely, if the cold entropy 
blobs were much smaller and were still to contribute signifi¬ 
cantly to the total energy flux, i.e., if they were much faster, 
this might not be compatible with an upward directed total en¬ 
ergy transport. A possible alternative might be suggested by 
the compressible simulations of Cattaneo et al. (1991), where 
the downward-directed kinetic energy flux in the downdrafts 
was found to balance the upward directed enthalpy flux in the 
downdrafts. This would imply that convection would trans¬ 
port energy only in the upwellings. This result has however 
not been confirmed in realistic surface simulations (Stein et 
al. 1992/ 

Equation (23) shows how fs decreases with depth, but its 
actual value and that of Vd remains undetermined. In the 
following, we propose a quantitative prescription for Vd by 
estimating its value within the top few hundred kilometers. In 
those top layers, we expect Vd to reach its maximum value, 
Vgiax, ]3e a certain fraction of V — Vad- In 

Appendix E we show that Vq®'’'/(V — Vad)max ~ 3. 

In deeper layers, where the local value of V — Vad has 
become small, Vd should scale with = f (AE)^, and 
therefore, using Equation (23), it should be proportional to 
fs{z) cx in addition to an Ma“^ factor; see Equa¬ 

tion (21). Thus, we have 

Vd = f! (32) 

^ Not to be confused with the parameter amix of Section 4.2 below. 
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where we have used for fs the expression 

fs = fsoip/p*)~^ (for p > p*). (33) 

Here, fso is a prefactor determining the strength of the result¬ 
ing Deardorff flux and p* is the density at the point at which 

V — Vad is equal to its maximum value (which is just below 
the photosphere). The new exponent ( = ( — A( takes the 
scaling of Mach number with density into account, i.e., 

Maocp-^S (34) 

where A( (> 0) will be computed in Section 4.3. 

3.4. Kinetic energy flux 

The importance of the kinetic energy flux has been stressed 
for some time as a property of compressible stratified convec¬ 
tion (Hurlburt et al. 1984; Cattaneo et al. 1991). This flux 
is related to the asymmetry of up- and downflows (Stein & 
Nordlund 1989) and is non-vanishing when / ^ 1/2; see 
Section 3.3. It is neglected in standard MLT, as has been dis¬ 
cussed by Arnett et al. (2015). Just like the Deardorff flux, 
the kinetic energy flux is also a non-gradient flux, but it is al¬ 
ways directed downward for / < 1/2 and must therefore be 
overcome by the enthalpy flux so that energy can still be trans¬ 
ported outwards. However, since the lowest order correlations 
in J^kin are triple correlations, there is no r approximation 
treatment analogous to that of the Deardorff flux. However, 
by comparing Fkin = —/)kinpu?ins with the enthalpy flux in 
Equation (20), it is possible to define a corresponding nabla 
term via Fkin = -^^cpT (TredW?n,s/^p)Vkin- This yields 

^kin — 3 ^kin'^rmsT/p/TpedCpT', (35) 

which is obtained analogously to Vd in Equation (21). The 
prefactor is here defined with a positive sign, so the total non- 
radiative flux caused by the turbulence is proportional to V — 

^ad ^kin- 

Obviously, if Vkin is strictly proportional to V — Vad, the 
addition of the Vkin term does not modify standard MLT, pro¬ 
vided that Vkin depends just on the local value of Urms, which 
in turn is, again, assumed to depend just on the local value of 

V — Vad- This is different for Vd, which depends on the en¬ 
tropy deficit produced by cooling near the surface and in this 
way on the value of V — Vad at the position where p = p,. 
The Vd term is therefore truly nonlocal. In the following, we 
combine kinetic and enthalpy fluxes into a total contribution 
Tconv = Tenth + Tkin which arises from convection. 

4. MODIFIED MIXING LENGTH MODEL 

4.1. Stratiflcation and flux balance 

To construct an equilibrium model, we begin by consider¬ 
ing first the case without convection, so the flux F is carried 
by radiation alone. Hydrostatic and thermal equilibrium then 
imply AP/Az = — pp and KAT/Az = —F, or, alternatively 
for the logarithmic gradients, 

dlnP/dz = -pp/P, (36) 

A\nT/Az = -F/{KT). (37) 

The double-logarithmic temperature gradient is obtained by 
dividing the two equations through each other, i.e., 

dlnT FP FcpWi^d 

V =-= = (38) 

dlnP KTpg Kg ' 


where we have used the perfect gas equation of state in the 
form P/Tp = cp — cy = cp (1 — I/7) = cp Vad- If the 
energy is no longer carried by radiation alone, V cannot be 
computed from Equation (38), but we have to invoke a suit¬ 
able theory of convection. In standard MLT, one obtains V as 
a solution of a cubic equation (Vitense 1953; Kippenhahn & 
Weigert 1990). In the following, we consider a modification 
that accounts for the possibility of a Deardorff flux. 

Elux balance implies that the sum of the radiative, enthalpy, 
and kinetic energy fluxes equals the total flux, i.e., 

-Ttot = Frad + ^conv ■ (39) 

where Fconv = T’enth + Tkin is the flux that arises from con¬ 
vection. In MLT, it is customary to express these fluxes in 
terms of nablas. Eor a given double-logarithmic temperature 
gradient V, the radiative flux is evidently 

Frad = V, (40) 

CpVad 

where V characterizes the actual temperature gradient. We 
can also define a hypothetical radiative temperature gradient 
Vrad that would result if all the energy were carried by radia¬ 
tion, so we can write 

Ptot = -§-Vrad, (41) 

CpVad 

which follows from Equation (38). We note in passing that the 
ratio between Ptot and the radiative flux carried by the adia¬ 
batic temperature gradient, Kg/cp, is known in laboratory 
and theoretical studies of convection as the Nusselt number 
(Hurlburt et al. 1984), whose local value is thus equal to 

NU = Vrad/Vad- (42) 

Einally, as explained in Section 3, we have from Equation (20) 

Penth = Po (V — Vad + Vd) , (43) 

with Po = |pcpPrredUrins/-Hp, and Pkin being essentially 
proportional to Penth; cf. Equations (30) and (31). Elux equi¬ 
librium then implies 

Vrad = V + e(V — Vad + VD), (44) 

where e = PoCpVad(l - (/kin//'enth)/-ff5- However, the 
expression for Fq involves the still unknown values of Wrms 
and Tred, which will be discussed in Section 4.2. 

4.2. Relaxation time and mixing length 

In convection the relevant time scale is the turnover time, 
which we write as t = l/itrmsfcf- We argue that kf should be 
estimated via the separation between the entropy rain struc¬ 
tures and not their thickness. This becomes important in the 
picture in which the downdraft threads of the entropy rain 
merge with neighboring ones to form a tree-like structure, 
as seen in the surface simulations (Stein & Nordlund 1989; 
Spruit 1997), leading therefore to different scalings of the 
two length scales (separation and thickness of structures) with 
depth. It is then not obvious which of these scales are more 
relevant for determining the kf that is relevant for mixing. In 
view of the uncertainty regarding the choice of kf, as well as 
for comparison with standard MET, we consider models with 
a fixed value of kf as well as the more conventional case in 
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Fig. 2.— Illustrative flow structures (upper row) and corresponding horizontal power spectra (lower row) associated with the three combinations of ^ and /3 
considered in this paper. Light shades coiTespond to large logai'ithmic power, which is seen to extend over lai'ge values of k for /3 = 0 (Cases I and II) and is 
confined to progressively smaller k at larger depths when /3 = 1 (Case III). 


which kfHp fv const. To capture the various cases in one 
expression, we assume in the following 

h = , (45) 

where /? = 0 corresponds to k{ = kiQ with a fixed value fcfo of 
the wavenumber of the entropy rain and /3 = 1 corresponds to 
kfHp = ttmix, where we have allowed for the possibility of 
a free mixing length parameter, amix. as is commonly done 
in standard MLT. It is not to be confused with the parameter 
omlt that was introduced in Section 3.3 and Appendix D. 
Negative values of /3 would correspond to shrinking scales, 
but will not be considered here. Likewise, non-integer values 
of [3 are conceivable, but will also not be considered here. 

Returning now to the discussion of the relaxation time r 
in the beginning of this subsection, instead of associating 
it with the turnover time we will allow for the 

possibility of an additional dilution factor (p^z) and write 
T = {u-cmskf4>)~^■ Here, 4>{z) increases with depth in a 
similar fashion as the scale height, so we assume, in anal¬ 
ogy with our treatment of kf in Equation (45), the expres¬ 
sion (p = (fcfoHp/ctmix)^”^- This dilution factor only en¬ 
ters in expressions involving the turbulent diffusivity, such as 
in Equation (17), and hence terms involving r or T^ed- For 
/3 = 1, the case (3 = 1 corresponds to the usual MLT concept, 
while ^ = 0 corresponds to a value of p that increases with 
depth. 


In Eigure 2 we present illustrative flow structures as well 
as their depth-dependent horizontal power spectra associated 
with the three combinations of (3 and (3 considered here. At 
the top of the domain, the size and separation of flow struc¬ 
tures are the same in all three cases. They remain constant 
with depth in the case (3 = (3 = 0 (Case I), while for (3 = 0 
and (3=1 (Case II), the separation increases with depth, 
but the thickness is still constant. Einally, for (3 = (3 = 1 
(Case III), both thickness and separation increase with depth. 
The Ailing factor decreases with depth in the case (3 = 0 
and /3 = 1, which could potentially make the kinetic energy 
flux divergent with depth, while for both (3 = (3 = 0 and 
(3 = (3 = 1, the Ailing factor is independent of height. 

With these preparations in place, we can write Fq in Equa¬ 
tion (43) in the form 

Fo = (a/37Vad)pUrmsc2a“fx(fcfoT?p)-(i-'^\ (46) 

where 

fj = Tj-edtirms^f^/^ — (ttrms F jp') (47) 

quantifies the radiative heat exchange between convective el¬ 
ements and the surroundings. The t term defined in Equa¬ 
tion (13), has a maximum at £kf = which typically oc¬ 
curs near the surface (e.g. Barekat & Brandenburg 2014). 

In Equation (46), the local value of the turbulent rms veloc- 
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ity always depends on the actual flux transported, and there¬ 
fore it must also depend on V — Vad + Vd- On dimensional 
grounds, since Fgnth is proportional to we have 

Wrms = Co (V — Vad + Vd)^^^ , (48) 

SO Fenth and Fkin are proportional to (V — Vad + Vd)^'^^- 
Standard mixing length arguments can be used to show that 
the prefactor in Equation (48) is a fraction of Cg. As shown in 
Appendix F, we have 

co/cs = V o'aMLT/37 \ (49) 

where /?' = (/? + /3)/2. It is then clear that Fkin scales with 
Hp like \ This is not the case for Fenth, 

however, because of the factor (fcfoFp)^“^ which, together 
with the (fcfoFp)^“^ factor, gives the scaling proportional to 
\ where /I" = (/3' + P)/2. 

4.3. Equation for the superadiabatic gradient 

Now that we know Fq, we can solve the equation for the su¬ 
peradiabatic gradient. This leads to an equation that is similar 
to the cubic equation for V, which is familiar from standard 
MET (Kippenhahn & Weigert 1990), 

Vrad = V + e* (V — Vad + Vd)^ , (50) 

where ^ = 3/2, and = (1 0kin//^enth)(centh F Ckin) 
has contributions from the enthalpy and kinetic energy fluxes. 
These expressions are similar to e in Equation (44), except 
that Eenth is evaluated with cq in place of Mrms, i-C-, 

Centh = (cr/37)3/2 /XQ) 

(51) 

where x = IXjpcp'vs, the radiative diffusivity. This shows that 
e* is essentially a Peclet number based on Cg. The contribution 
from the kinetic energy flux is given by 

t^kin/t^enth ~ /^kinC-MLT^ad ^ ■ (52) 

We note in passing that e* is also related to the Rayleigh 
number, which is commonly defined in laboratory and nu¬ 
merical studies of convection; see Appendix G. Furthermore, 
because of convection and the resulting bulk mixing, S is 
now approximately constant, and therefore, unlike in the non- 
convecting reference solution with K = const (Section 2), 
K can no longer be constant, but it reaches a minimum at the 
point where k is maximum, which turns out to be at a depth of 
about 1 Mm in the convection models presented below. Since 
e* is inversely proportional to K, it reaches a maximum at 
that depth and falls off both toward the top and the bottom of 
the convection zone. 

An essential difference between Equation (50) and the usual 
one in MET is the presence of Vd arising from the Dear- 
dorff flux. Within the usual MET, where Vd = 0, one finds 
that V is slightly above Vad, but now it might instead be 
slightly above Vad — Vd- There are indeed two possibili¬ 
ties for convecting solutions (Fconv > 0), one corresponding 
to a Schwarzschild-stable solution, 

Vad - Vd < V < Vad (stable), (53) 

and one that is Schwarzschild unstable, 

Vad - Vd < Vad < V (unstable). (54) 
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Fig. 3.— Solution of Equation (50) for 5 = 3/2 and V^ad = 10® (solid 
black), compared with the approximation (56) for q = (dashed blue) 
and g = 1 (dash-dotted red), and ^ = 1 (thin orange). The dotted line gives 
an additional unphysical solution of Equation (50) for ^ = 3/2. The limiting 
cases V = Vad ^nd V^ad ^re shown as thin horizontal green lines. 


Which of the two possibilities is attained depends on the value 
of Vd and also on details of the solution. As will be discussed 
in Section 6 below, entropy rain convection may actually still 
be Schwarzschild unstable without exciting giant cell con¬ 
vection if small-scale turbulent viscosity and diffusivity are 
strong enough so that the local turbulent Rayleigh number for 
the deeper layers is subcritical. 

In this connection, we note that in standard MET, one in¬ 
cludes the effects of radiative cooling of the convective ele¬ 
ments in a different manner than here. Instead of V — Vad, 
the effective buoyancy force is written as V — V', where V' 
always lies between V and Vad (Vitense 1953). Thus, one 
has Vad < V' < V, which resembles Equation (54) with a 
negative value of Vd- 

To understand the nature of the solutions of Equation (50), 
it is instructive to treat ^ as an adjustable parameter. For given 
values of Vrad and Vad = Vad — Vd, the case ^ = 1 yields 

\ Vrad + C*Vad 

V £*) =-—-- (55) 

1 + e* 

It shows that V Vrad for e* —> 0 (stable surface layers) and 

V —Vad for e* ^ 1 (deeper layers). Next, to discuss the 
general case ^ ^ 1, we define AV = V — Vad and AVrad = 
Vrad — Vad- For AVrad < 0 We have AV = AVrad, while 
for AVrad > 0 and AV 1, a useful approximation is 

AV AV,'A / (gAV;A-' + , (56) 

where 9 = ft agrees with Equation (55) in the special 
case ^ = 1, where AV = AVrad/(l + e*)- In Figure 3 we 
plot V versus e* for Vrad = 10®- The approximation yields 

V > Vrad for e* < 10“®, which is unphysical. This can be 
mitigated by choosing q = 1', see Figure 3. 

For the relevant case of large values of e*, we have 

AVf« (Vrad/e*)'A. (57) 

This relation is useful because, even though both Vrad and e* 
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TABLE 2 

Comparison of Cav . 2C = — 2AC, and their difference for 

VARIOUS COMBINATIONS OF /3 AND 0 USING C = 1. 



/3 

P 

p' 

P" 

Cav 

2C 

Cav - 

Case I 

0 

0 

0 

0 

4/9 

2/9 

2/9 

Ca.se II 

0 

I 

1/2 

3/4 

10/9 

2/9 

8/9 


1 

0 

1/2 

1/4 

2/3 

2/3 

0 

Case III 

1 

1 

1 

1 

4/3 

2/3 

2/3 


depend on K, their ratio does not and is given by 

Vrad 3Vad ^tot/p (1+2/3") 

-=- 7 ^ oc --ocT 

e* (T cog/k{ pci 

(58) 

where /3" = (/3+3/3)/4. Thus, AV oc T gives the 

scaling of A V for the bulk of the convection zone. Therefore, 

looking at Equation (48), we find Ma oc T ™ with 

m = 1 -/3'+ (1/2 +/3")/e (59) 

For ^ = 3/2, using the relation 3/3' — 2/3" = /3, we find 
that m = (4 — ^)/3 is independent of the value of /3. Thus, 
we have m = 4/3 with /3 = 0 and m = 1 with /3 = 1. 
For isentropic stratification, this implies for the Mach number, 
given by Equation (34), the following scaling: A(/ = 8/9 with 
/3 = 0 and AC = 2/3 with /3 = 1. 

4.4. Relative importance of the Deardorff term 

In Section 3.3 we have considered the depth dependence of 
the Deardorff term via Equations (32) and (33). However, for 
Vd to be important at increasing depths, it must exceed the 
sufiadiabatic gradient Vad — V, because otherwise it would 
not be possible for the Deardorff term to make V — Vad + 
Vd positive. From Equations (57) and (58) we see that AV 
depends on T in a power-law fashion. We are particularly 
interested in the conditions under which AV falls off faster 
than Vd, because that would ensure that the Vd term remains 
important even at larger depths. 

For C = 3/2, and since the convection zone is nearly isen- 
tropically stratified (T ex we have 

AVoep-'^^^ with Cav = ^(1 + 2^"). (60) 

In Table 2 we compare for various combinations of /3 and /3 
the exponents for AV and Vd, i-e., 2C = 2C and 2AC, where 
AC = 2m/3. In all cases we have chosen C = 1, i-e., we 
allow for moderately non-ideal (radiative) effects relative to 
the ideal case with C = 0.8. We see that the difference is 
positive and non-vanishing in all cases, except for /3 = 1 and 
^ = 0. In the following, we present solutions for all of the 
remaining three cases. Before doing this, let us recapitulate 
what led to the threefold dominance of (3 over /3. For better 
illustration, we summarize in Table 3 the various relationships 
that led to the scaling of AV with kf^Hp. 

In the first expression for fenth, /3 enters because it char¬ 
acterizes the relation between the buoyancy force propor¬ 
tional to ST/T and advection proportional to see Ap¬ 

pendix D. For thin threads, we expect the relevant kf to be 
large, i.e., /3 = 0 (Cases I and II in Figure 2). Next, in the 
second expression for FCnth, we have used a mean-field ex¬ 
pression to relate Fenth to AV via a turbulent diffusivity pro¬ 
portional to ft! Urms/kfcj), where, as argued above, the 


TABLE 3 

Illustration of scaling relationships with kfoHp . 


^enth IX ^rms {fioRpf ^ 

Fnth ^ rirms ASJ / [kfoHpf-P 

OC /{kfoHpf-P' 

Fenth OC (AV)3/2/(fcfoFp)2(l-/5") 
Fkin OC (AV)3/2/(fcfoFp)3(l-P') 


buoyancy force 
mean-field expression 

0 ' = { 0 + 0)/2 

0" = {0'+~0)/2 = {0+W)/A 

kinetic energy flux 


dilution factor </> has entered. It is this expression that is clos¬ 
est to conventional MET, because here we expect /3 = 1. 

The remaining two relationships in Table 3 explain why the 
/3 term appears three times more dominantly than the /3 term. 
By equating the first two expressions for Fenth in Table 3, we 
find first of all the relation between itrms and AV, where /3 
and /3 contribute with equal shares through /3' = (/3 -f $)/2. 
However, to see the scaling of AV, we need to go back to the 
second expression for Fenth, because it changes only weakly 
with depth. Now, /3' and /3 contribute with equal shares, and 
this means that /3 has now become three times more dominant 
than /9 through the expression /3" = (/3 -f 3^)/4. Thus, for 
/3 = 0 and /3 = 1, AV shows nearly the standard scaling with 
Hp. Furthermore, looking again at the first expression for 
Fenth in Table 3, we see that only (3 enters, so the scaling of 
Urms is fully characterized by that of small blobs with negative 
buoyancy. 


5. NUMERICAL SOLUTIONS 

In this section we present numerical solutions to demon¬ 
strate the effect of the Vd term on the resulting stratifica¬ 
tion. At the end of this section, we also compare with the 
non-convecting reference solution mentioned in the introduc¬ 
tion. We should emphasize that, although we use solar param¬ 
eters, our models cannot represent the Sun, because ioniza¬ 
tion effects have been ignored (we take /i = 0.6 for the mean 
molecular weight in Cp — Cv = TZ/p, where TZ is the uni¬ 
versal gas constant). A rather simple opacity law of the form 
of Equation (6) with kq = lO^cm^g"^, po = 10“^gcm“^, 
Tq = 13, 000 K, and a = 0.5, b = 18, for the exponents in the 
power-law expression for kh- has been used. We also neglect 
the departure from plane-parallel geometry, so our model can 
only give qualitative indications. 

The system of two differential equations (36) and (37) de¬ 
couples by using In P as the independent variable. We thus 
integrate 

dlnT/dlnF = V, (61) 

using Equations (32), (33), and (50) to compute V. As 
an initial condition we use Ftop = Teff/2^/^ at a suffi¬ 
ciently low pressure (here Ftop = lO^dyncm”^), so as 
to capture the initially isothermal part of the atmosphere; 
see, e.g., Bohm-Vitense (1958) or Mihalas (1978). Here, 
Teff = (Ftot/cse)^^'^ is the effective temperature. We ap¬ 
proximate the Itrms term in Equation (47) by using the value 
from the previous step. The Deardorff term is characterized 
by the assumed value of fso and the value of — A(/ 

(Section 3.3), where C = 1 A^ is a function of (3 and (3 
(Section 4.3). Since we integrate from the top downward, no 
prior knowledge of p* is needed, because the Deardorff term 
is invoked only after V — Vad has reached its peak value. 
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Fig. 4.— Profiles of Sjcp, V — Vad, and Urms for fso = 0 (Vd = 0), 
as well as fso = 0.2, 0.3, and 0.5 for ^ ^ = 0 (Case I) with ^ = 1/9. 

The location of the surface (r = 1) is indicated by vertical dash-dotted lines 
and geometric depths below the surface are indicated in the middle panel, 
stalling with tick mai'ks at 100, 200, and 500 km, and continuing with 1, 
2, and 5 Mm, etc. The inset in the middle panel shows V — Vad over a 
nari'ower range as a function of —z. 

Geometrical and optical depths are obtained respectively as 

—z = J{P/pg)d\nP and = J (62) 

where the integration of r starts at InPtop and that of — z at 
the position where r = 1, which is referred to as the surface. 
The factor kP/ g, which is the same as Hp, is retained be¬ 
cause of the similarity with that in the expression for r. The 
z coordinate is used in some of our plots. 

In the following, we present solutions for the three combi¬ 
nations of [3 and [3 sketched in Figure 2. We recall that only 
Cases I and II (/3 = 0 with /3 = 0 or 1) correspond to small 
length scales in the deeper layers (see also the power spectra 
in Figure 2) and are therefore of interest when trying to rec¬ 
oncile the non-detection of convective motions by Hanasoge 
et al. (2012, 2016) at the theoretically expected levels. We did 
already emphasize that Case II with /? = 0 and /3 = 1 is likely 
to lead to large kinetic energy fluxes. However, based on the 
results presented below, it turns out that in the case /3 = /? = 0 
(Case I), which was favored by these two requirements (small 
length scale and non-divergent kinetic energy flux), the Dear- 
dorff term is unlikely to have a significant effect, because for 
/3 = 0, the gradient term in the enthalpy flux becomes rather 
inefficient and must therefore be compensated for by a cor¬ 
respondingly larger superadiabatic gradient, and thus, only a 





Fig. 5.— Same as Figure 4, but for /3 = 0, /3 = 1 (Case II) with (^ = 1/9, 
and /jjo = 0, 0.1, and 0.2. 

rather large Deardorff flux (/so > 0.5) can make the resulting 
stratification sufficiently subadiabatic. This is the case shown 
in Figure 4, where we present profiles of iS/cp, V — Vad, and 
Urms for /so = 0 (no Deardorff term), as well as /^o = 0.2, 
0.3, and 0.5. The Urms profiles are basically the same for 

all values of fso and fall off like P . This is expected, 

because Urms c)c and m = 4/3; see Equation (59), 

— 1/2 

where we have used Cs oc T . Thus, for isentropic stratifi¬ 
cation we have Urms oc 2m)/5 _ p 1/3^ 

Next, we show in Figures 5 and 6 Cases II and III with (3 = 
1 and /so = 0, 0.1, and 0.2. For both f3 = 0 and f3 = 1, the 
Deardorff flux now has a stronger effect and is able to make 
the deeper parts of the domain Schwarzschild stable even for 
fso = 0.1. Those layers would then be Schwarzschild stable 
and no longer a source of giant cells. Again, the profiles of 

Urms are similar regardless of the value of fso, but fall off 

_ 

more slowly for 13 = 1 (Urms oc P ), compared to /3 = 0 
(Urms oc P ^^^). This agrees with our theory, because for 
m = 1 we have Urms oc P*'^ 2m)/5 _ p 1/5^ 

Given that the stratification is stable in the deeper parts, we 
can calculate the Brunt-Vaisala frequency of buoyancy oscil¬ 
lations, Nby, which is given by iVgy = —(V — Va,d)g/Hp- 
At intermediate and larger depths of the convection zone, we 
have g/Hp ^ 10“^ and Vad — V ^ 10“^, so the pe¬ 
riod of buoyancy oscillations would be of the order of days. 
This is comparable to or less than the turnover time r. Indeed, 
one finds that tNbv ~ \/Vd/AV — 1 exceeds unity in the 
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Fig. 7.— Comparison of 5/cp (top) and Vrad/Vad (bottom) between the 
non-convective radiative reference solution (V = V^ad) and standai'd con¬ 
vective solutions (Vd =0) with /3 = /3 = 0 (red, dashed) and 0 = 0 = 1 
(blue, dotted). In both panels, the location of the r = 1 surface is indicated 
by vertical dash-dotted lines and geometric depths below the surface are in¬ 
dicated for the non-convective solution, starting with tick marks at 50, 100, 
and 200 km, etc. The location of the a priori unstable layer (V^ad > Vad) 
is marked in the lower panel by a small gray strip. For the other solutions, 
the depths are different; see those in Figure 6 for the solution with fso = 0, 
which is the same as here for 0 = 0 = 1. 


Fig. 6.— Same as Figure 4, but lor 0 = 0 = 1 (Case III) ij = 1/3, and 
fsO = O-lj and 0.2. 

deeper parts, which is a consequence of Vd falling off with 
p with a smaller power than AV, as discussed in Section 4.4. 
Nevertheless, the resulting decrease in S/cp with depth re¬ 
mains always small (10“^ or less) compared with the value 
of AS'o/cp produced at the surface (« Mrms/®MLTVadCs ~ 
0.01). Thus, based on this, the descending low-entropy blobs 
would reach the bottom of the convection zone before their 
negative buoyancy is neutralized by the decreasing average 
entropy. In other words, they will never perform any actual os¬ 
cillations before reaching the bottom of the convection zone, 
i.e. where T/onv = 0. 

It may be interesting to note that the depth of the con¬ 
vection zone increases slightly with increasing values of /so- 
Of course, our model is idealized and represents the Sun at 
best only approximately. Furthermore, as emphasized in Sec¬ 
tion 1, the depth of the convection zone is well determined 
seismically, and this should be reproduced by a solar model 
with realistic atomic physics and appropriately chosen ad¬ 
justable parameters. However, it is known from the work of 
Serenelli et al. (2009) that a slight expansion of the solar con¬ 
vection zone would actually be required to compensate for the 
shrinking that follows from the downward revision of solar 
abundances (Asplund et al. 2004) which is based on three- 
dimensional convective atmosphere simulations, compared to 
previous analysis based on one-dimensional semi-empirical 
models (Grevesse & Sauval 1998). 

Finally, we compare in Figure 7 the standard convective so¬ 
lution (/3 = /? = 1; same as the case with /so = 0 in Figure 6) 
with the non-convective radiative reference solution. Not sur¬ 
prisingly, owing to the absence of convection, the same flux 


can now only be transported with a greatly enhanced nega¬ 
tive (unstable) entropy gradient near the surface. However, 
this layer is now extremely thin (1.15 Mm) and the peak in 
Vrad/Vad IS about 100 km below the t = 1 surface. Note 
also that its peak value (~ 10®) is below that in the presence 
of convection, where it reaches a maximum of ~ 4 x 10® at 
a depth of « 1 Mm. As shown in Equation (42), this value 
can be interpreted as the local Nusselt number. Note also that 
the result for /3 = /3 = 0 (when Vd is weak) is rather similar 
to that for /3 = ^ = 1; see the dashed and dotted lines in 
Figure 7. 

Our calculations have demonstrated that for fi = 1, regard¬ 
less of the value of j3, bulk mixing changes the non-convecting 
reference state to a nearly isentropic one. However, whether 
the mean entropy gradient is slightly stably or slightly unsta¬ 
bly stratified depends on the presence of the Deardorff flux. 
In the model with ^ = 0, however, bulk mixing is rather inef¬ 
ficient and the stratification would be Schwarzschild unstable 
unless an unrealistically large Deardorff flux is invoked. 

At the end of the introduction, we discussed that the entropy 
rain itself might create an unstable stratification. Let us now 
return to this question with more detailed estimates. This will 
be done in the following section. 

6. ALTERNATIVE CONSIDERATIONS 

Assuming that the upflows are perfectly isentropic. Spruit 
(1997) argues that the low-entropy material from the top with 
its decreasing entropy filling factor fs in deeper layers neces¬ 
sarily leads to a negative mean entropy gradient. Specifically, 
using Equation (24), one obtains 

= ( 63 , 

dz cp din P 5 cp 
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Fig. 8.— Dependence of U^g/uj-^s on depths — 2 : for Case I with /so = 
0.5 (solid red line), as well as Cases II and III with /so = 0.1 (dashed blue 
and dotted green lines, respectively). 

where we have used d In /g/dlnp = 2/3 and din p/d In P = 
3/5 for an isentropic layer with 7 = 5 / 3 . Thus, the stratifica¬ 
tion would be Schwarzschild unstable. This is also borne out 
by the solar simulations with realistic physics, although the 
computational domains are sufficiently shallow so that the ra¬ 
diative flux is still small in the deeper parts and usually even 
neglected altogether. Toward the bottom of the convection 
zone, however, radiation becomes progressively more impor¬ 
tant and the mean entropy gradient in the upflows may no 
longer be vanishing. 

To estimate the mean entropy gradient in the upflows, we 
may balance the steady state entropy advection with the neg¬ 
ative radiative flux divergence, i.e. 

pTU^dS^/dz fa-{dFra,d/dz)^. (64) 

The sign of dPrad/dz is negative and thus compatible with 
a positive dSf jdz in the upflows, but it would only be large 
enough near the bottom of the convection zone. Higher up, 
dPrad/d^; becomes smaller and eventually unimportant. On 
the other hand, C/-|-/Mrnis can be rather small (see Table 1). 
Furthermore, our considerations neglect the fact that the gas in 
the upflows expands, so only a fraction of the gas can ascend 
before it begins to occupy the available surface area. There¬ 
fore, the rest of the gas would have to remain stagnant and 
continue to heat up. In reality, of course, there would be con¬ 
tinuous entrainment, resulting in a finite, but still low, effec¬ 
tive upward velocity. 

We now use selected solutions obtained in Section 5 to es¬ 
timate the effective fractional upward velocity for radiative 
heating/cooling to dominate over advection, defined as 

C/Iff/Urms = (-dPrad/dz) / (Urms P TdS'/dz) . (65) 

Here we have used TdSjdz = p AV/Vad- Figure 8 shows 
this quantity for Cases I-III, which are shown in Figures 4-6. 
It turns out that for Case I with /^o = 0.5, the effective frac¬ 
tional upward velocity is around 0.01, while for Case III with 
/so = 0.1 it can be even larger. The value 0.01 is compatible 

with the C/Vurms values in Table 1 if we assume / = 0.01 

(so C/Vurms = 0.1) and [/gfj/[/^ = 0.1. For Case II with 
/sO = 0.1, the effective fractional upward velocity is around 
10“^, which may be unrealistically small. 

Based on the considerations above, we may conclude that 
the case for a Schwarzschild-stable entropy gradient depends 


on model assumptions for the implementation of the Dear- 
dorff flux that are potentially in conflict with the mean en¬ 
tropy gradient expected from the two-stream model. Whether 
or not there is really a potential conflict depends not only on 
the model parameters, but also on Equation (64) itself. This 
is because the T factor in the entropy equation is outside 
the derivative, making it impossible to derive the total flux 
balance assumed in Equation (39); see also Rogachevskii & 
Kleeorin (2015) for related details. 

In view of these complications, it is worthwhile to discuss 
alternative ways of avoiding giant cell convection. Eor this 
purpose, we have to address the global problem of using some 
coarse-grained form of effective mean-field equations. When 
considering the equations of mean-field hydrodynamics, in 
which the small-scale enthalpy and momentum fluxes are pa¬ 
rameterized in terms of negative mean entropy and mean ve¬ 
locity gradients, respectively, one finds solar differential rota¬ 
tion as a result of non-diffusive contributions to the Reynolds 
stress, but in certain parameter regimes a new instability was 
found to develop (Tuominen & Rudiger 1989; Rudiger 1989; 
Rudiger & Spahn 1992). This instability was later identified 
as one that is analogous to Rayleigh-Benard convection, but 
now for an already convecting mean state (Tuominen et al. 
1994). This instability would lead to giant cell convection. 

The existence of giant cell convection is under debate, but 
assuming that it does not exist in the Sun, one might either 
hypothesize that the mean-field equations are too simplified 
or that mean-field convection could be suppressed by suffi¬ 
ciently strong turbulent viscosity and turbulent thermal diffu- 
sivity coefficients. These turbulent coefficients define a turbu¬ 
lent Rayleigh number for a layer of thickness d, 


Ra. = ^ (f) 


( 66 ) 


This number would then have to be still below the critical 
value for convection. Equation (66) differs from the usual 
one defined through Equation (Gl) in Appendix G in that, 
first, V ^ vt and x ^ Xt have been substituted, and sec¬ 
ond, the superadiabatic gradient of the non-convecting refer¬ 
ence solution is now replaced by the actual one. This might 
be a plausible alternative to explaining the absence of giant 
cell convection if the idea of turning the stratification from 
Schwarzschild unstable to Schwarzschild stable through the 
Vd term were to turn out untenable. 

We recapitulate that in this alternate explanation, the strat¬ 
ification is Schwarzschild unstable, i.e. Rat > 0, correspond¬ 
ing to Equation (54), but Rat < Ra“'* is still below a certain 
critical value, so it would be stable by a turbulent version of 
the Rayleigh-Benard criterion. Tuominen et al. (1994) found 
Ra“'* ft! 300 for a vanishing rotation rate. However, they also 
estimated that Rat > Ra“'* for plausible solar parameters, 
so the direct adoption of this idea would be problematic, too. 
However, one might speculate that a more accurate treatment 
could lead to stability when allowing, for example, for spa¬ 
tial nonlocality of the turbulent transport (Brandenburg et al. 
2008; Rheinhardt & Brandenburg 2012). 


7. CONCLUSIONS 


In the present work we have suggested that the enthalpy flux 
in stellar mixing length models should contain an extra non¬ 
local contribution so that the enthalpy flux is no longer pro¬ 
portional to the local superadiabatic gradient, V — Vad, but 
to V — Vad + Vd, where Vd is a new nonlocal contribution 
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that was first identified by Deardorff (1966) in the meteoro¬ 
logical context. The significance of this term lies in the fact 
that it provides an alternative to the usual local entropy gradi¬ 
ent term and can transport enthalpy flux outwards—even in a 
slightly stably stratified layer. 

We have presented a modified formulation of stellar MLT 
that includes the Vd term, in addition to a Vkin term result¬ 
ing from the kinetic energy flux. The formalism and the fi¬ 
nal results are similar to those of conventional MLT in that 
one also arrives here at a cubic equation for V, but the term 
V — Vad is now replaced by V — Vad + Vd — Vkin- This 
new formulation implies that convection can carry a finite flux 
while V — Vad is still negative and therefore the stratification 
is Schwarzschild stable, i.e.. Equation (53) is obeyed. Conse¬ 
quently, if confirmed, no large length scales are being excited. 

The present formulation allows for different treatments of 
the length scales governing buoyant elements on the one hand 
and the time and length scales associated with mixing on the 
other. When both are independent of depth (/3 = /3 = 0), 
mixing becomes inefficient at larger depths. Thus, to carry a 
certain fraction of the enthalpy flux, the superadiabatic gradi¬ 
ent needs to be larger than otherwise, making it harder for the 
Deardorff term to revert the sign of V — Vad- On the other 
hand, for a tree-like hierarchy of many downdrafts merging 
into fewer thin ones at greater depths (/3 = 0, /3 = 1), 
the increasing length scale associated with increasing separa¬ 
tion enhances vertical mixing, making the stratification nearly 
isentropic without the Deardorff term, and slightly subadia- 
batic with a weak Deardorff term. In that case, however, the 
Ailing factor of the downflows decreases with depth as p~‘', 
which may imply an unrealistically large downward kinetic 
energy flux. This leaves us with the standard flow topology 
(/3 = /3 = 1), where both size and separation of structures 
increase with depth. The Deardorff flux can still cause the 
stratification to have a subadiabatic gradient, so no giant cell 
convection would be excited locally, but the flow structures 
would be large and should be helioseismically detectable, as 
has been found by Greer et al. (2015) using ring-diagram local 
helioseismology. 

It would be useful to explore the thermodynamic aspects of 
the present model more thoroughly and to connect with re¬ 
lated approaches. An example is the work by Rempel (2004), 
who studied a semianalytic overshoot model that was driven 
nonlocally by downdraft plumes, similar to what was sug¬ 
gested by Spruit (1997). Rempel also finds an extended suba¬ 
diabatic layer in large parts of the model. Furthermore, there 
are similarities to the nonlocal mixing length model by Xiong 
& Deng (2001), who, again. And an extended deeper layer 
that is subadiabatically stratified in the deeper parts of their 
model. In this connection we emphasize the main difference 


between nonlocal turbulence owing to the Deardorff flux and 
usual overshoot; in the latter case the enthalpy flux would 
go inward as a consequence of the reversed entropy gradient, 
while in the present model the Deardorff flux goes outward. 

Future work could proceed along two separate paths. On 
the one hand, one must establish the detailed physics leading 
to the Vd term using models with reduced opacity, in which 
reliable DNS are still possible, i.e., no SGS terms are added 
and the primitive equations are solved as stated, without in¬ 
voking Equation (1). On the other hand, one could study suit¬ 
ably parameterized large eddy simulations that either include 
a nonlocal Deardorff term of the form given by Equation (32), 
as discussed in Section 3, or that explicitly release entropy 
rain at the surface such that the resulting stratification is still 
slightly stable. This would be particularly useful in global 
simulations that would otherwise have no entropy rain. 

As stimulating as the results of Hanasoge et al. (2012) are, 
they do require further scrutiny and call for the resolution of 
the existing conflicts with other helioseismic studies such as 
those of Greer et al. (2015). Alternatively, global helioseismic 
techniques for detecting giant cell convection (Lavely & Ritz- 
woller 1993; Chatterjee & Antia 2009; Woodard 2014) can 
provide another independent way of detecting deep larger- 
scale flows (Woodard 2016). Realistic simulations should 
eventually agree with helioseismic results of flows in deeper 
layers of the Sun, but at the moment it is still unclear whether 
a subgrid scale treatment as in Equation (1) adequately cap¬ 
tures the small-scale flows that can be responsible for the 
Deardorff flux and whether they would in principle be able 
to predict the subtle departures from superadiabatic stratifica¬ 
tion on subthermal time scales. 
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APPENDIX 

A. POLYTROPIC STRATIFICATION FROM KRAMERS OPACITY 

We show here that, for the non-convecting reference solution, using the Kramers-like opacity law of Equation (2), but not the 
combined opacity law of Equation (6), we have K const in the deeper, optically thick layers. Dividing Equation (3) by the 
equation for hydrostatic equilibrium, AP/dz = —pg, we have 

dr _ frad _ -Frad {p/PoT _ fiad jP/P qT 
dP “ Kpg ~ Kopog {T/nr-b - Kopog (T/To)3+“-^’ 


where Kq = IGassT^ / {3kqPq) is a constant and P/Pq = {p/ pq){T/Tq) is the ideal gas equation with a suitably defined 
constant Pq = (cp — Cv)poTo. Here, po and Tq are reference values that were defined in Equation (2). Equation (Al) can be 
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integrated to give 

(r/ro)"+“-'’ = {n+ l)v(°^(P/Po)'+“ + (Ttop/To)4+“-^ (A2) 

where = Fra-dPo/(KoToPog), which is defined analogously to the Vrad without superscript (0) in Equations (38) and (41), 
and Ttop is an integration constant that is specified such that T —>■ Ttop as P —>■ 0. Note also that 4 + a — 6 = (n + 1)(1 + a), 
where n was defined in Equation (5) as the polytropic index, so the ratio of 4 + a — 6 to 1 + a is just n + 1, which enters in front 

of the term in Equation (A2). Since K oc T^-^j cx T^+^-^j we have K —>■ const = Kq for T ^ Ttop. 


B. DERIVATION OF EQUATIONS (10) AND (II) 

To obtain Equations (11) and (12), which are used to derive the Deardorff flux term in the r approximation, we start with the 
equations for specific entropy and velocity in the form (see, e.g., Barekat & Brandenburg 2014) 


BU 


— V • Frad, 

-VP + pg, 


(Bl) 

(B2) 


where D/Df = d/dt + (7 • V is the advective derivative, Frad is the radiative flux, and viscosity has been omitted. Subtracting 
those equations from their averaged ones, we obtain the following set of equations 


ds jj 1 

dt ^ dxj dxj pT 


V • bFrad 


+ A4, 


dui — dUj 

dt ^ dxj dxj 


= -V + 9i{p/lP - s/cp) +Afui, 
P 


(B3) 

(B4) 


where we have used 5p/'p = p/^P — s/cp, which can be obtained from Equation (8) and the perfect gas law by linearization. 
Pressure fluctuations will again be neglected and (V • i5Frad)/pT will be replaced by —s/tcooI, as explained in Section 3. 
Assuming U = 0 and omitting the nonlinear terms A4 and Afui we arrive at Equations (11) and (12). 

C. FILLING FACTOR FOR A DESCENDING HILL VORTEX 

The solution for a Hill vortex with radius an and propagation velocity uh is given by a stream function rp in spherical coordi¬ 
nates (r, 6, (j)) as (e.g. Moffatt &. Moore 1978) 


1 f —3 uh(1 — VO (for r < or), 

4 \ +2'Uh(1 - Oh/''’^) ^ (for ^ > “h). 


(Cl) 


We apply it in Cartesian coordinates as the initial condition for the mass flux as pw = V x ('!'</>), where cf) = {—y/vu^xjw, 0) 
is the unit vector in the toroidal direction of the vortex, using w'^ = x^ + and = vu'^ + for cylindrical radius w and 
spherical radius r. We adopt here an isothermal equation of state, i.e., there is no buoyancy force in this problem. For non- 
isothermal calculations, but in two dimensions, we refer to the work of Rast (1998). Density and pressure fall off exponentially 
with height and both the scale height and the sound speed are independent of height. No analytic solution exists in that case, so 
the Hill vortex solution is at best approximate. We consider a domain of size L x L x AL with —Lj^ < x, y/Hp < L/2 and 
—7T/2 < z < Ll2, where L = hPlp. We choose or = O.SFp and ur = 0.2cs. The viscosity is z/ = 5 x 10“®, so the Reynolds 
number is appaY^jv = 2000. We use the PENCIL CODE^ with a resolution of 1152 x 1152 x 4608 meshpoints. 

Figure 9 shows snapshots zoomed into the vortex as it traverses about five scale heights. The filling factor, which is proportional 
to the radius squared, decreases with depth and is found to scale with the surrounding density like Equation (23) with (j = 0.8; 
see the last panel of Figure 9. 


D. MLT RELATION BETWEEN Urms AND Srms 

To find the relation between itrms and the rms values of temperature or entropy fluctuations, it is customary in standard MLT 
to approximate the steady state momentum equation, m • Vtt « —g ST/T, by = omlt g ST/T, where omlt ~ 1/8 is a 

commonly adopted geometric factor (e.g.. Spruit 1974). This leads to 

Srms/Cp = ST/T = Mrms^f/®MLTff = T/p/oMLT, (Dl) 

where we have used Cg = ^gHp, and thus (srms/cp) Ma“^ = 'ykfHp/aMLT- Furthermore, using VadCpT = gHp, we have 
Tsrms = tt?ms^fffp/(®MLTVad), which is used to derive Equation (31). 

^ https://github.com/pencil-code 


17 



Fig. 9.— Velocity vectors superimposed on a color scale representation of the vorticity uj at times 0, 10, and 20 in units of {Hp /as well as the resulting 
filling factor versus density. The negative slope is f = 0.8. Note that the frame of view changes as the vortex descends and shrinks. 


E. ESTIMATE FOR THE SURFACE VALUE OF Vd 

The purpose of this section is to show that Vd is a certain fraction of V — Vad in the top layers, as stated in Equation (32). To 
have an estimate for we multiply Equation (11) by s and average, so we get 

1 0s‘^ _ _ 

- s'^/Tcoo\+Ts. (El) 

As for fenth in Equations (14) and (15), we have a triple correlation term, which is here Ts = —su ■ Vs. Again, we adopt the r 
approximation and replace 7^ by a damping term of the form Ts = —s^jr which, together with the TcooI term, combines to give 
Tred as the relevant time scale. Assuming a statistically steady state, ds'^/dt = 0, we derive the following expression for s^: 

s^ = — Tred (in the surface layers). (E2) 

This shows that fluctuations of specific entropy are produced when there is an outward flux (UTs > 0) and a locally negative 
(unstable) mean entropy gradient; see also Garaud et al. (2010) for a similar derivation. Inserting this into Equation (18), and 
using Equation (9), we obtain 

Fb = g (Fenth • VS/cp) (in the surface layers). (E3) 

Here, both g and 'VS point downward, so Fp, points upward and we can write Fp, = AFenth, where A = (Tred/'nr)^(V — Vad) 
is a coefficient that itself is proportional to the superadiabatic gradient, and r® = (Hp/g)^/^ is the free-fall time (after which a 
fluid parcel at rest has reached a depth of Hp/2). Since Fenth = Fq + Fd, this implies Fq = (1 — A)Fenth and therefore A < 1 
in the highly unstable layer at the top, where both Fq and Fgnth are positive. We expect Fd to be largest just a few hundred 
kilometers below the photosphere, so A should be maximum in the upper parts. To calculate the fraction Vd/(V — Vad), we use 
the fact that Fenth = Fq + Fd, together with the part of Equation (20) that relates to Fd, to write 

Fd = \'pCpT (rredUrms/-^F')VD = (Tred/Tff)^ (V — Vad)Fenth- (E4) 

Starting with Equation (31) and expressing Tred in terms of a using Equation (47), we And 

Vr7(V - Vad)max = Sui^enthVad a-^ihoEp)-^^-^^ ^ 3, (E5) 

where we have used /3 = 1 (appropriate to the near-surface layers), ignored radiative cooling (cr = 1), and assumed amix = 1.6 
(Stix 2002), as well as (penth ~ 4 that was found in simulations of Brandenburg et al. (2005), as discussed in Section 3.3. 

We emphasize that these equations only characterize the initiation of entropy rain. They cannot be used to compute the 
Deardorff flux in the deeper layers where we have instead invoked Spruit’s concept of nonlocal convection in the form of threads. 
The associated with those threads is likely to result from a part of the triple correlation term Ts = -su-'Vs, which gives 
rise to a negative divergence of the flux of of the form H = us^ on the rhs of Equation (El). This flux (which should not 
be confused with an energy flux) should point downward (s^ is largest when < 0) and be strongest in the upper layers, so 
V H < 0, leading to a positive contribution from — V • iF on the rhs of Equation (El). This may explain the associated with 
the Vd term. 

To estimate the ratio Fd/Fq in the deeper layers, we use Equation (El) in the steady state with I/tcooI —S' 0 and T = 
— V • H Ki (2C — A()s'^Urms/sHp and obtain 

0 = —UjsVjS + (2( — A()s'^u^^s/lHp (in deeper layers). (E6) 

Multiplying by irredWrms (P^)^, using Equations (9) and (17), and expanding the fraction by g/cp, we have 

0 = Fenth Fg + ip7ms(2C “ ^C) (rreds'^9PT^ CpT/{'jgHp). (E7) 
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This implies that Fq < 0 in the deeper layers. The second term in the parentheses is Fd (with a plus sign, because p > 0 is a 
scalar here); see Equation (18). Furthermore, we use Cpr/( 7 pi?p) = 1 /( 7 — 1) and = T'enth/</i’enth; see Equation (31). 

The Fenth terms on both sides cancel, so we have 

Fb/\Fg\ = 3 ( 7 - l)/[</.enth(2C-AC)] ~ 0.3... 0.5, (E8) 

where 2( — AC = 2( + AC =10/9 for Cases I and II, and 14/9 for Case III (Section 4.3), and (penth = 4 has been assumed. 


F. DERIVATION OF EXPRESSION FOR co 

To find the coefficient cg given by Equation (49), we use Equation (Dl) in the form 

Cp6T = U^^s/CfiFp/VadOMLT- 

Inserting this into Equation (7), and using Equation (45), yields 

Fenth = pUrms “mix(^fO-^T')^~^/^adaMLT- 

Equating this expression with Equation (43), using Equation (46), we can derive the desired expression for rtrms in the form 


'l£ ^ _ 


flMLT 


3 ^ (V - Vad + Vd) , 


3 7 (fcfoFp)2-2/5 

where /3' = (/3 + /3)/2 and Cg = jgHp have been used. We thus find the coefficient cg as stated in Equation (49). 


(FI) 

(F2) 

(F3) 


G. RAYLEIGH NUMBER 

The purpose of this section is to show that e,, as defined in Equation (51), is related to the Rayleigh number Ra. In laboratory 
and numerical studies of convection, it is customary to define Ra as (e.g. Kapyla et al. 2009) 


^ gd^ gj^ / V.ad-Vad y°"~^°"^,_. gj^ / V.ad-Vad \ 

V ^ V Hp ) vx \ Hp ) 


which is usually evaluated in the middle of the layer at z = 2 *. In the Sun, however, the maximum value is more relevant. The 
superscript “non-conv” indicates that the entropy gradient is taken for the non-convecting reference state, d is the thickness of the 
layer and p is the viscosity. Introducing the Prandtl number Pr = p/x and using the definition of Nu in Equation (42), we have 


PrRa _ gd'^Vad 
Nu - 1 “ x'^Hp ■ 

On the other hand, using Equation (51), we find 

2 _ g^QMLT i.cl/xa)'^ _ O'^QMLT Fpff _ O'^QMLT (Fp/d)"^ ffd^^Vad 

“ 2773 (fcfgFp)4(l-/3") - 27(fcfgFp)4(l-F') x^ “ 27Vad(A:fgFp)4(l-/3") x^Hp 


(G2) 


(G3) 


and therefore 

2 _ cr3aMLT/27Vad PrRa 

(fcfod)4(fcfoFp)-4/3" Nu-r 


(G4) 


which shows that Ra is proportional to ej. 
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